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We construct a theoretical framework to describe the evolution of heavy flavors produced in 
relativistic heavy-ion collisions. The in-medium energy loss of heavy quarks is described using our 
modified Langevin equation that incorporates both quasi-elastic scatterings and the medium-induced 
gluon radiation. The space-time profiles of the fireball are described by a (2-|-l)-dimensional hydro¬ 
dynamics simulation. A hybrid model of fragmentation and coalescence is utilized for heavy quark 
hadronization, after which the produced heavy mesons together with the soft hadrons produced from 
the bulk QGP are fed into the hadron cascade UrQMD model to simulate the subsequent hadronic 
interactions. We find that the medium-induced gluon radiation contributes significantly to heavy 
quark energy loss at high pt; heavy-light quark coalescence enhances heavy meson production at 
intermediate px; and scatterings inside the hadron gas further suppress the D meson Raa at large 
Pt and enhance its V 2 . Our calculations provide good descriptions of heavy meson suppression and 
elliptic flow observed at both the LHC and RHIC. 


I. INTRODUCTION 

The primary purpose of performing relativistic heavy- 
ion collisions at the Large Hadron Collider (LHC) and 
the Relativistic Heavy-Ion Collider (RHIC) is to study 
the properties of QCD matter under extreme conditions 
such as high temperatures and densities. It has now been 
well established that a new state of matter, known as the 
strongly interacting quark-gluon plasma (sQGP) Qi , is 
created in these energetic nuclear collisions. This highly 
excited state of matter is composed of color de-confined 
quarks and gluons, and displays properties of a nearly 
perfect fluid such as the strong collective flow observed 
for the produced hadrons [^-Q- Relativistic hydrody¬ 
namic models successfully describe the space-time evolu¬ 
tion of the strongly coupled QGP fireballs nn , from 
which it is found that the value of the shear viscosity to 
entropy density ratio p/s of the produced QGP is small. 

Apart from studying soft hadrons emitted from the 
QGP directly, an alternative way to study the transport 
properties of the QGP is through the investigation of 
the modification to the properties of energetic partons 
that travel through the produced hot and dense medium. 
One of the promising candidates is heavy quarks. Ow¬ 
ing to their large masses, the thermal production of 
heavy quarks from the QGP fireball is significantly sup¬ 
pressed, thus the majority of them are produced dur¬ 
ing the primordial stage of the collision via hard scat¬ 
terings. These heavy quarks then propagate through the 
medium, and can probe the whole evolution history of the 
QGP. Over the past decade, experimental observations 
at both the LHC and RHIC have revealed many interest¬ 
ing and sometimes surprising observations of heavy flavor 
hadrons and their decay electrons, such as the small val¬ 
ues of their nuclear modification factors Raa and the 
large values of their elliptic flow coefficients which are 


almost comparable to those of light hadrons (l^-[l9l| . This 
seems contradictory to the earlier expectation from the 
mass hierarchy of parton energy loss and still remains a 
challenge for us to fully understand. 

Various transport models have been constructed to 
study the heavy quark motion inside dense nuclear mat¬ 
ter, such as the par ton cascade model based on the Boltz¬ 
mann equation (20l - l23l | and the linearized Boltzmann ap¬ 
proach coupled to a hydrodynamic background [l^, [l^ ■ 
In the limit of small momentum transfer, the multiple 
scatterings of heavy quarks inside a thermalized medium 
can be treated as Brownian motion, and the Boltzmann 
equation for quasi-elastic scatterings is then reduced to 
the Fokker-Plank equation which can then be stochasti¬ 
cally realized by the Lan gevi n equation. Many Langevin- 
based transport models [2fl-[^ have been developed to 
study the collisional energy loss of heavy quarks and have 
been shown to be successful in describing experimental 
data in the low transverse momentum px region where 
the phase space for the medium-induced gluon radiation 
is restricted by the lar ge m asses of heavy quarks, i.e., 
the “dead cone effect” [351 Is^ . However, LHC exper¬ 
iments now enable us to observe heavy meson spectra 
up to 30 GeV. At such high px, even heavy quarks be¬ 
come ultra-relativistic and therefore the radiative energy 
loss should no longer be neglected. In our previous work 
(3^ [3^ , the classical Langevin equation is modified such 
that quasi-elastic scattering and medium-induced gluon 
radiation can be incorporated simultaneously. In this 
study we will continue utilizing this improved Langevin 
approach for the in-medium evolution of heavy quarks. 

A dedicated description of the heavy quark energy loss 
inside the QGP is crucial for solving the “heavy flavor 
puzzle”, but has yet to be accomplished. It has been 
pointed out that details in the hadronization process may 
have a strong impact on the observed heavy meson spec- 
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tra [ 2 ^, The influence of hadronic interactions 

after the QGP decays on heavy meson observables has 
also been explored in Refs. [13, and has been shown 
to be non-negligible. In this work, we will develop a 
hybrid model of fragmentation plus coalescence to de¬ 
scribe the heavy quark hadronization process. The mo¬ 
mentum dependence of the relative probability between 
fragmentation and coalescence will be calculated accord¬ 
ing to the Wigner functions in an instantaneous coales¬ 
cence model. This coalescence model was first proposed 
for the production of light hadrons out of QGP fireballs 
(idl - li^ , and then applied to the production of heavy fla¬ 
vor hadrons in nuclear collisions and recently to 

partonic jet hadronization as well. This model does 
not require the thermalization of the recombining par- 
tons and is easily extensible to simultaneously include 
various meson and baryon species, allowing for the nor¬ 
malization of the total coalescence probability over all 
possible hadronization channels. Based on our previous 
study [S^ , we will further develop this hybrid hadroniza¬ 
tion model such that it is applicable to arbitrary local 
flow velocities of the QGP background. In addition, the 
rescattering of heavy mesons inside a hadron gas after 
hadronization will also be incorporated in this work by 
utilizing the Ultra-relativistic Quantum Molecular Dy¬ 
namic Model (UrQMD) , and the effect of hadronic 
interactions on the observed heavy meson spectra will be 
investigated in detail. 

The paper is organized as follows. In Sec. im we will 
present how the classical Langevin equation is modified 
to simultaneously incorporate collisional and radiative 
energy loss of heavy quarks and how the simulation of 
the heavy quark evolution in a dynamic QGP medium is 
implemented. In Sec. nil we develop a hybrid model of 
fragmentation and coalescence to describe the hadroniza¬ 
tion of heavy quarks. With that hadronization model, 
we present numerical results of heavy meson suppression 
and anisotropic flow and compare to experimental data 
at the LHC and RHIC. In Sec. IVl we will discuss how 
the hadronic rescattering of heavy mesons is simulated 
within the UrQMD model and its effect on the observed 
heavy meson Raa and vi . We will summarize and discuss 
future developments in Sec. 


II. HEAVY QUARK ENERGY LOSS IN QGP 
MATTER 

A. A modified Langevin equation 

During their propagation through a thermalized QCD 
matter, heavy quarks lose energy via both quasi-elastic 
scatterings with light patrons in the medium and gluon 
radiation induced by multiple scatterings. In this work, 
we utilize the following modified Langevin equation [s^ 
that simultaneously incorporates these two processes to 
describe the time evolution of energy and momentum of 


heavy quarks while they traverse the QGP matter: 

= -m{v)v + i + fa- ( 1 ) 

In Eq. ([TJ, the first two terms on the right-hand side 
are inherited from the classical Langevin equation and 
represent the drag force and the thermal random force 
experienced by a heavy quark while it diffuses inside a 
thermal medium due to multiple scatterings. For a mini¬ 
mal model, we assume the thermal force ^ is independent 
of the heavy quark momentum and satisfies the correla¬ 
tion relation of a white noise it')) = — t'), 

in which k denotes the momentum diffusion coefficient of 
heavy quarks and is related to the spatial diffusion coeffi¬ 
cient D via D = T/[M77 d(0)] = 2T^/k if the fluctuation- 
dissipation theorem rjjjij)) = k/{2TE) is respected. 

Apart from the above two forces resulting from quasi¬ 
elastic scatterings, an additional term fg = —dpg/dt 
is introduced into Eq. to describe the recoil force 
exerted on heavy quarks while experiencing medium- 
induced gluon radiation, where Pg is the momentum of 
the radiated gluon. The probability of gluon radiation 
during the time interval [f, t + At] is determined based 
on the average number of radiated gluons in this At: 

Pi.ad(t, At) = (Vg(t, At)) =^t J (2) 

As long as At is chosen sufficiently small, (Ag(t, At)) 
is less than I and can be interpreted as a probabil¬ 
ity. In this study, the gluon distribution function in Eq. 
([2]) is adopted from the higher-twist calculation for the 
medium-induced gluon radiation - the distribution func¬ 
tion of gluons radiated from a massless parton is calcu¬ 
lated in Refs. [H, and its modification due to the 
mass effect of a heavy quark is introduced by Ref. : 

dNg _ 2asP{x)q . 2 ( t - tA ( k\ ^ 
dxdk^dt Tvk'^ \ 2Tf j \k‘^+x'^M'^J ’ 

in which x is the fractional energy taken by the emit¬ 
ted gluon from the heavy quark, and k± is the trans¬ 
verse momentum of the gluon, as is the strong cou¬ 
pling constant, P{x) is the gluon splitting function and 
Tf is the formation time of the gluon defined as Tf = 
2Exil — x)/{k‘]_ x^M^) with E and M being the en¬ 

ergy and mass of heavy quarks. Note that the multi¬ 
plicative term at the end of Eq. @ is known as the 
“dead cone factor”, signifying the mass dependence of 
the radiative energy loss of hard parton. In Eq. ©, 
q is the gluon transport coefficient and may be related 
to the above mentioned quark diffusion coefficient k via 
q = 2 kCa/Cf- Therefore, in our calculations there is 
only one free parameter in the modified Langevin equa¬ 
tion [Eq. ® ]. To obtain the best description of heavy 
flavor observables at the LHC and the RHIC, as will be 
shown in Sec. m and Sec. lYl the spatial diffusion co¬ 
efficient of heavy quark D{2'kT) is chosen around 5^6, 
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which is equivalent to qjT^ around 9.4 ^ 11.3 for glu¬ 
ons (or 4.2 ~ 5.0 for quarks), consistent with the value 
extracted by the JET Collaboration via fitting the exper- 
imental data using various jet energy loss models . 

When simulating the radiative energy loss of heavy 
quarks, a lower cut-off of radiated gluon energy wq = ttT 
is imposed to take into account the balance between 
gluon emission and absorption processes. Below ojq, the 
gluon radiation is disabled and the evolution of heavy 
quarks with low energies is entirely controlled by quasi¬ 
elastic scatterings. For this reason, x G [ttT/E, 1] is used 
when calculating the gluon radiation probability in Eq. 
©• With this treatment, we have verified in our pre¬ 
vious work that the thermal equilibration of heavy 
quarks can be approached after sufficiently long evolu¬ 
tion time although the exact fluctuation-dissipation rela¬ 
tion may not be guaranteed due to the lack of the gluon 
absorption process. A more detailed discussion of this 
approach and possible improvements towards a more rig¬ 
orous treatment of detailed balance between gluon emis¬ 
sion and absorption were discussed Reference, [s^. And 
alternative approaches of including the radiative energy 
loss into the Langevin framework can be found in Refs. 

fsEIH. 


B. Heavy quark evolution in a realistic medium 

To study the heavy flavor spectra produced in real¬ 
istic heavy-ion collisions, we couple the above modified 
Langevin equation to an expanding QGP medium that 
is simulated with a (2-|-l)-dimensional viscous hydrody¬ 
namic model developed in Refs. [^, [Tl|, [l^. In this work, 
we utilize the code version and parameter values provided 
by Ref. 0. The hydrodynamic simulation generates the 
space-time evolution of the local temperature and flow 
velocity profiles of the QGP fireball created in relativis¬ 
tic nuclear collisions. For every time step of the Langevin 
evolution, we first boost each heavy quark into the local 
rest frame of the fluid cell through which it propagates. 
In the rest frame of fluid cell, the energy and momentum 
of a given heavy quark are updated using Eq. dm before 
it is boosted back to the global center of mass frame. 

The hydrodynamical evolution of the bulk matter is 
initialized with either the Monte Carlo (MC) Glauber 
or the Kharzeev-Levin-Nardi (KLN) parametrization of 
the Color Glass Condensate (CGC) model for its en¬ 
tropy density distribution. To best describe the spec¬ 
tra of soft hadrons emitted from the QGP fireballs, for 
both the RHIC and the LHC environments, the starting 
time of the QGP evolution has been set as tq = 0.6 fm/c 
and the shear-viscosity-to-entropy-density ratio (rf/s) has 
been tuned as 0.08 when the Glauber initial condition is 
used and 0.20 when KLN is used. In this work, a smooth 
initial condition is utilized for the bulk matter. Possible 
effects of the the initial state fluctuation on heavy fla¬ 
vor observables have been discussed in our earlier study 
[ 5 ^ . For heavy quarks, we use the MC Glauber model 




FIG. 1: (Color online) The initial heavy flavor spectra from 
the leading-order pQCD calculation with and without the nu¬ 
clear shadowing effect (EPS09), (a) for the LHC and (b) for 
the RHIC experiments. 


to initialize their production positions and the leading- 
order perturbative QCD (pQCD) approach to cal¬ 
culate their initial momentum space distribution. We 
have included the pair production process {gg —>■ QQ and 
qq —>■ QQ) and the flavor excitation process {gQ —?► gQ 
and gQ -G gQ) in calculating the initial px spectra of 
heavy quarks. The gluon splitting process {g —>■ QQ) 
has been recently discussed in Ref. and will be in¬ 
vestigated in a follow-up study. These pQCD calcula¬ 
tions are at the partonic level. To calculate the cross 
sections of heavy quark production in nuclear collisions, 
we adopt CTEQ for the parton distribution functions 
and include the nuclear shadowing/anti-shadowing effect 
in heavy-ion collisions using the EPS09 parametrization 


In Fig. [U we show the px spectra of initial heavy 
quarks at both LHC and RHIC energies, for proton- 
proton collisions and binary collision number scaled 
nucleus-nucleus collisions. The influence of the nuclear 
shadowing/anti-shadowing effect in the initial state on 
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heavy quark spectra can be clearly observed in the fig¬ 
ures: it reduces the production rate of charm quarks at 
low pt but slight enhances it at high px; the effect is 
stronger at the LHC energy than at the RHIC. For bot¬ 
tom quarks, the production of the low px bottom quarks 
is decreased at the LHC energy but slightly enhanced at 
the RHIC when initial state effects are included. Such 
effects will have significant impact on the nuclear modi¬ 
fication factor Raa of heavy mesons observed in the fi¬ 
nal state as will be shown later in Sec. ttm The cal¬ 
culated spectra are used to sample the initial px distri¬ 
butions of heavy quarks. Their initial rapidity distribu¬ 
tions are taken to be uniform around the central region 
(-1 < ?7 < 1 ). 

In simulating the evolution of heavy quarks, they are 
assumed to stream freely first from their production ver¬ 
tices in hard collisions to tq = 0.6 fm/c, the initial time 
at which the hydrodynamical evolution commences. The 
possible energy loss in the pre-equilibrium stage has been 
neglected which is expected to give only a small contri¬ 
bution to the final state spectra, given its short period of 
time compared to the much longer evolution of the QGP 
fireball. 

With the above setup, we can investigate how heavy 
quarks evolve inside QGP matter and lose their energy. 
In Fig. O we calculate the average energy loss of charm 
and bottom quarks as a function of their initial energy 
after they traverse a realistic QGP medium created by 
central Pb-Pb collisions at the LHC energy. The contri¬ 
butions from different energy loss mechanisms are com¬ 
pared. As shown by the figures, for both charm and 
bottom quarks, quasi-elastic scatterings dominate their 
energy loss while their initial energies are small, however, 
medium-induced gluon radiation dominates in the high 
energy regimes. The crossing point is around 7 GeV for 
charm quarks, and increases to 18 GeV for bottom quarks 
due to the greater suppression of gluon radiation by its 
larger masses. These results indicate that collisional en¬ 
ergy loss alone may provide reasonable descriptions for 
the heavy flavor observables in the low px region as those 
measured at RHIC, but will become insufficient when we 
extend to higher px such as those reached by the LHC 
experiments. 


III. HEAVY FLAVOR HADRONIZATION 

In the previous section, we studied the initial produc¬ 
tion of heavy quarks in heavy-ion collisions and their en¬ 
ergy loss inside a QGP medium. Around the critical 
temperature Tc = 165 MeV, both the bulk matter of the 
QGP fireball and heavy quark should hadronize into color 
neutral bound states. For the bulk matter, we utilize the 
numerical tool “iSS” [i^ based on the Cooper-Frye for¬ 
mula to obtain soft hadrons from the hydrodynamic 
medium. For heavy quarks, we follow our previous work 
and develop a hybrid model of fragmentation and 
coalescence to describe their hadronization. After heavy 




FIG. 2: (Color online) Comparison of energy loss between 
different mechanisms: (a) for charm quark and (b) for bottom 
quark. 

mesons are obtained from heavy quarks, we may directly 
compare their suppression and collective flow coefficients 
with experimental data from both the LHC and RHIC. 


A. A hybrid model of fragmentation and 
coalescence 

Between the two typical in-medium hadronization pro¬ 
cesses, fragmentation and heavy-light quark coalescence, 
of heavy quarks into heavy flavor hadrons, the former 
dominates the high momentum regimes while the latter 
becomes important at low momenta. The momentum 
dependence of the relative probability between these two 
mechanisms can be determined by the Wigner function 
in the instantaneous coalescence model [^. With the 
knowledge of this probability, spectra of heavy mesons 
formed from the heavy-light quark coalescence can be 
directly calculated within the coalescence model, while 
those from the fragmentation process can be obtained 
from the Pythia simulation [66|. In our previous study 
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(3^ . a hybrid model of fragmentation and coalescence 
was established for the heavy flavor hadronization. In 
this work, we will further develop this hadronization 
model so that the effect of the local flow of an expanding 
medium on hadronization can be conveniently taken into 
account. 

In the instantaneous coalescence model, the momen¬ 
tum spectra of produced mesons and baryons are given 
as follows. 


The above procedure can be straightforwardly general¬ 
ized to a three-body system for baryon formation by com¬ 
bining two quarks first and then combining their center 
of mass with the third quark: 




with qi and q 2 as the relative momenta defined in the 
rest frame of the produced baryon 


dNM 

(PpM 

dNe 

d^PB 


= / d'^pid'^p2 


-3 dNi dN2 


-I 


d^pi d?p2'' ^ 


j3 dNi dN2 dN^ yy 


= / d'^pid-^p2d'^p3- 


d?pi d^p2 d^p3 
xS{pM -Pl-P2 -Ps)- 


fM{'Pl^P2)5{pM - Pi -P 2 ) 
{PUP2,P3) 


(4) 


dNi/d^Pi denotes the momentum distribution of the i- 
th valence quark in the produced hadron. The spectra of 
heavy quarks can be directly obtained after they traverse 
the QGP fireball within our modified Langevin evolution. 
Light quarks are assumed thermal in the local rest frame 
of the expanding medium: 


dN, ^ ggV 

d^P ^y/p'^+m^/T _|_ 2 ’ 


(5) 


in which pq = 6 is the statistic factor that takes into 
account spin and color degeneracy of quark, and for sim¬ 
plicity a uniform distribution in the position space is as¬ 
sumed inside a volume V. In Eq. o, one key ingredi¬ 
ent of the coalescence model is the Wigner function 
which denotes the probability for the two or three quarks 
to combine. For a two-body system, the Wigner function 
can be written as 

fM{r,^ = gM j d^r'e~^‘>''' + ( 6 ) 


in which qm denotes the degrees of freedom (spin and 
color) of the formed meson and the variables f and q are 
the relative position and momentum of the two particles 
defined in the two-body center-of-mass frame, i.e., the 
rest frame of the produced meson: 


77'cm^cm ZT" cm,^cm 

^2 Pi ^1 P 2 


Ecm 




(7) 


Note that the heavy and light quarks are first boosted 
into their center-of-mass frame in which their coalescence 
probability is then calculated. In Eq. (pM represents 
the meson wavefunction, which is approximated by the 
ground state wavefunction of a simple harmonic oscil¬ 
lator: exp[—r^/(2(T^)]/(7r(7^)^/'^. Here, the width a is 
related to the angular frequency of the oscillator uj via 
a = 1/y/jmjj, with p = mim 2 /{mi -|- m 2 ) being the re¬ 
duced mass of the two-body system. With these setups, 
we may average over the position space of Eq. © and 
obtain the momentum space Wigner function of the pro¬ 
duced meson: 


/m(9^) = 9m 


(20Fcr)3 g 2^2 


( 8 ) 


77’cm.;^cm 77 ’cm.^cm 

^2 Pi - P 2 


9i = 




Er 


77’cm/^cm I ^cm\ /ir'cm i 77'cm\^ 

- _ -^3 (Pi +P2 )- (El + E2 )P' 
92 = - 


■cm 




Er 


rcm 

^3 


, ( 10 ) 


and (7i = 1/y/piUj as the width parameter with pi = 
77117712 /{"mi + m 2 ) and p 2 = {mi -f 7712 ) 7773/(7771 -|- 7772 -I- 
7773 ). In the calculations, the thermal mass is taken as 
300 MeV for u and d quarks and 475 MeV for s quarks. 
Heavy quarks, on the other hand, are not required to be 
thermal, and their masses are taken as 1.27 GeV for c and 
4.19 GeV for b quarks. Gontribution from thermal gluons 
is also incorporated in this coalescence model: they are 
split into light quark pairs first and then combine with 
heavy quarks to form heavy flavor hadrons. 

We use Eqs.® and m to evaluate the momentum 
dependence of heavy-light quark coalescence probabili¬ 
ties at the critical temperature Tc, as shown in Fig. |31 
In principle, the oscillator frequency w in these Wigner 
functions can be calculated from the charge radius of 
the hadrons and should depend on the hadron species. 
Here for a minimal model, we adopt an average value 
0.215 GeV for all c-hadrons and 0.102 GeV for 6-hadrons. 
These two parameters are obtained by requiring the co¬ 
alescence probability through all possible hadronization 
channels to be unity for a zero momentum heavy quark 
in a static medium at Tc since it is not sufficiently ener¬ 
getic to fragment as can be seen in Fig. [3] In our 
calculations, all major hadron channels are incorporated, 
including the ground states and the first excited states 
of D/B mesons, Ag, Sg, Sg and Hg. 

After the to parameters are evaluated in a static 
medium according to the above normalization procedure, 
the Wigner functions are determined. For heavy quarks 
in an expanding medium, we adopt an effective temper¬ 
ature method to calculate the effective tempera¬ 

ture of a fluid cell with a non-zero velocity due to a blue 
shift effect as follows: 



g-^g>s/^eff J- 


V f d^v __ 

^IE ± 1 ’ 


( 11 ) 


in which u is the 4-velocity of the fluid cell. This ef¬ 
fective temperature is then utilized in the thermal 
distributions of light partons and the coalescence proba¬ 
bility inside a moving fluid cell is calculated according to 
Eqs. m and ([2]). If the obtained value of the coalescence 
probability is greater than unity at low momenta, it is 
taken as unity. 


V 
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PHQ(GeV) 

(a) 



FIG. 3: (Color online) The momentum dependence of the 
coalescence probabilities at different flow velocities: (a) for 
charm quark and (b) for bottom quark. 




(b) 

FIG. 4: (Golor online) Gomparison between the contributions 
from different hadronization mechanisms to the (a) D and (b) 
B meson spectra (normalized to one heavy quark). 


In Fig. [3l we show our calculations of the coales¬ 
cence probabilities for both charm and bottom quarks as 
functions of their momenta either through heavy meson 
channel alone {D/B meson), or to any possible hadrons 
(summing over all hadron channels under consideration). 
Three different values of the fluid flow velocity, which cor¬ 
respond to three different effective temperatures are com¬ 
pared. One can observe that the coalescence probability 
generally decreases with the increase of heavy quark mo¬ 
mentum, and a larger fluid velocity leads to a higher effec¬ 
tive temperature and therefore an enhanced coalescence 
probability. Furthermore, for the same momentum, bot¬ 
tom quarks have larger probability to coalesce with light 
quarks than charm quarks do, due to the larger mass 
(or smaller velocity) of the bottom quarks inside a QGP 
medium. 

In Fig. [3] we divide the hadronization of heavy quarks 
into three regimes: coalescence with light quarks to D 
or B mesons, coalescence to other hadron channels, and 
fragmentation. After its evolution through the QGP 
matter, if a charm or bottom quark is selected for coales¬ 


cence into a D or B meson, a light quark or anti-quark 
is generated according to thermal distribution at Tes in 
the local rest frame of the fluid cell, and then boosted 
to the lab frame to combine with the given heavy quark 
according to the probability governed by Eq. (|8l). If they 
do not combine, another light parton is generated until 
a meson is formed. On the other hand, if a heavy quark 
is selected to fragment based on the probability in Fig. 
[21 its fragmentation is implemented via Pythia in which 
the relative ratios between different hadron channels are 
properly calculated and normalized. 

Using this hybrid model of hadronization, we may com¬ 
pare the relative contributions from coalescence and frag¬ 
mentation mechanisms to heavy meson production in rel¬ 
ativistic nuclear collisions. As can be observed in Fig. |4l 
after charm and bottom quarks traverse a realistic QGP 
medium created in central Pb-Pb collisions at the LHC 
energy, their hadronization to D/B mesons are domi¬ 
nated by fragmentation at high px but is significantly 
enhanced by heavy-light quark coalescence at interme¬ 
diate pt- Since the coalescence mechanism combines a 
thermal parton and a heavy quark, the spectrum of D/B 
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mesons is shifted to the larger momentum regime com¬ 
pared to the original charm/bottom quark distribution. 
Therefore, its contribution to the production of heavy 
mesons at low px is not as significant as that at interme¬ 
diate pt- Furthermore, as already seen in Fig. [Sj due to 
the larger masses and thus smaller velocities of b quarks 
than c quarks, the coalescence mechanism dominates a 
wider px range for B meson production than for D me¬ 
son production. 


B. Heavy flavor suppression and collective flow 


With our modified Langevin equation for the in¬ 
medium evolution of open heavy quark and the above 
hybrid model of fragmentation and coalescence for heavy 
quark hadronization, we are able to calculate the suppres¬ 
sion and elliptic flow coefficients of heavy flavor hadrons 
and compare them with experimental data from the LHC 
and RHIC. Discussions on the additional variation of the 
heavy flavor observables due to the hadronic interactions 
after the QGP freezes out will be deferred to the next 
section. 

Because of the medium modification, heavy flavor 
hadrons produced in nucleus-nucleus collisions display 
different spectra from those produced in proton-proton 
collisions. The two most widely utilized quantities that 
characterize the medium effect are the nuclear modifica¬ 
tion factor Raa and the elliptic flow coefficient V 2 - 


Raa{pt) = 


1 dN^^/dpT 
N^i dNPP/dpT ’ 


n 2 (px) = (cos(2(^)) 


/ pI-pI 

\pI+pI 


( 12 ) 

(13) 


which describe the overall energy loss and the asymmetric 
Px modification of the probe particles respectively. 

In Fig. 5(a) we show our calculation of the D meson 
Raa in central Pb-Pb collisions at the LHC energy. The 
impact of the nuclear shadowing effect in heavy quark 
production in the initial state and the contribution of the 
coalescence mechanism to the D meson formation can be 
clearly observed in the figure. As shown in Fig. |5(a)l if 
other factors are fixed, the inclusion of the initial state 
shadowing effect would lead to a factor of 2 suppression 
of the D meson Raa at low px and a mild enhancement 
at high Px. This is consistent with the findings shown in 
Fig. 1(a) the production of charm quark is significantly 
suppressed at low px and slightly enhanced at high px 
in Pb-Pb collisions compared to that in proton-proton 
collisions. Therefore, a better understanding of the cold 
nuclear matter effect in the initial state is crucial for a 
more precise description of the heavy flavor suppression 
in nuclear collisions. From Fig. 5(a)[ we also observe 
that although the fragmentation mechanism alone is suf¬ 
ficient for describing the heavy quark hadronization at 
highpx (above 8 GeV), the coalescence of light and heavy 




FIG. 5: (Color online) The D meson (a) Raa and (b) V2 in 
2.76 TeV Pb-Pb collisions [H, [l^, compared between differ¬ 
ent hadronization mechanisms, and between with and without 
the nuclear shadowing effect. 


quarks becomes crucial in the low and intermediate re¬ 
gion: it converts low px heavy quarks into intermediate 
Px hadrons by combining the former with thermal par- 
tons from the QGP medium, and thus suppresses the D 
meson Raa near zero px but greatly enhances it in be¬ 
tween 2 and 5 GeV. With the incorporation of the nuclear 
shadowing effect in the initial state, a modified Langevin 
equation that includes both collisional and radiative en¬ 
ergy loss of heavy quarks inside the QGP matter, and 
a hybrid model of fragmentation and coalescence, our 
calculation provides a good description of the D meson 
Raa in central Pb-Pb collisions as measured by the AL- 
IGE Gollaboration. The spatial diffusion coefficient of 
heavy quark is determined as 5/(2ttT) by comparing our 
calculation to experimental data at high px, and will be 
utilized for all the following calculations in this section. 

Figure [F(b)| shows our results of the D meson V 2 in pe¬ 
ripheral Pb-Pb collisions at the LHC. The nuclear shad¬ 
owing effect is included for all the curves shown in this 
figure, but various hadronization mechanisms are com- 





































pared in more details. For the pure fragmentation pro¬ 
cess, the Wigner function /'^ in the coalescence model is 
set as a constant 0 in order to switch off all coalescence 
channels; to the contrary, /'^ is fixed at 1 for the pure co¬ 
alescence hadronization. One can observe that the pure 
coalescence limit leads to a much larger D meson V 2 than 
the pure fragmentation limit because the former mecha¬ 
nism brings the anisotropic flow of light quarks from the 
hydrodynamic background into the formation of heavy 
mesons. However, only a slight enhancement in the D 
meson V 2 at intermediate px is observed in our hybrid 
hadronization model compared to the pure fragmenta¬ 
tion process despite the large enhancement of its yield. 
This may result from the momentum dependence of the 
Wigner function in this instantaneous coalescence model 
that prefers combining partons with similar velocities. 
Other factors may also affect the final D meson V 2 such 
as the initial heavy quark spectra and the development 
of the radial flow in the hydrodynamic background. 

This may result from a combinational effect of the ini¬ 
tial heavy quark spectra, the momentum dependence of 
the Wigner function in this instantaneous coalescence 
model, and the development of the radial flow in the 
hydrodynamic background. 

In Fig. 0 we study the suppression and the elliptic 
flow of D mesons produced in the RHIC experiments. 
Although at the RHIC energy, the nuclear shadowing ef¬ 
fect for the low px heavy quark is not as significant that 
at the LHC energy, it still has a non-negligible impact 
on the D meson Raa as shown in Fig. 6(a) Since the 
current RHIC experiments concentrate on the relatively 
low px region, the introduction of heavy-light quark coa¬ 
lescence is even more crucial in the hadronization process 
than that for describing the LHC data. The coalescence 
mechanism results in a bump structure of the D meson 
i?AA around 1-2 GeV, which cannot be obtained with the 
pure fragmentation mechanism. In Fig. 6(b) we can see 


that the introduction of the coalescence mechanism helps 
increase D meson V 2 , similar to the findings in the LHC 
scenario. By including all the effects discussed above, our 
numerical results are consistent with the STAR data. 


One of the most interesting puzzles related to heavy 
flavor is the mass hierarchy of parton energy loss. In Fig. 
[71 we compare the suppression between D and B mesons 
due to different energy loss mechanisms, in which the 
mass hierarchy of heavy quark energy loss can be clearly 
observed for both quasi-elastic scattering and medium- 
induced gluon radiation. Due to their larger masses, 
bottom quarks lose significantly smaller amount of en¬ 
ergy than charm quark does at low px after they prop¬ 
agate through a realistic QGP medium and therefore B 
meson displays larger i?AA than D mesons. With our cur¬ 
rent model calculation, the mass effect on collisional en¬ 
ergy loss becomes negligible for the meson spectra above 
20 GeV. However, difference in radiative energy loss still 
remains up to 40 GeV. Apart from the mass hierarchy of 
the in-medium parton energy loss, there is also the mass 
dependence for heavy-light quark coalescence probability 



4 

Pt (GeV) 
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FIG. 6: (Color online) The D meson (a) Raa and (b) V2 
in 200 GeV An-An collisions [3, [13) compared between dif¬ 
ferent hadronization mechanisms, and with and without the 
nuclear shadowing effect. 


as shown in Fig. [3 Since it is easier for bottom quarks 
to combine with thermal partons from the medium back¬ 
ground than for charm quarks, the enhancement of the 
B meson Raa is more prominent than that of the D me¬ 
son Raa ; such enhancement also spreads over a wider px 
regime for B mesons. 

One possible direct verification of the mass hierarchy of 
parton energy loss is the comparison of the nuclear sup¬ 
pression of D mesons versus non-prompt J/tp as shown 
in Fig. [5] Here we show our calculations of the partic¬ 
ipant number dependence of the Raa for D mesons, B 
mesons and non-prompt J/ip. The decay from B meson 
to J/tp is implemented with Pythia. As has been men¬ 
tioned earlier, the only free parameter in our transport 
model is the spatial diffusion coefficient of heavy quarks 
which is fixed to be D = 5/(2TrT) by comparing high px 
D meson Raa in 2.76 TeV central Pb-Pb collisions to 
experimental data. One can see that with a single value 
for the transport coefficient, our calculation provides a 
good description of the participant number dependence 
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FIG. 7: (Color online) Heavy meson suppression in cen¬ 
tral Pb-Pb collisions, compared between different energy loss 
mechanisms and D and B mesons. 



FIG. 8: (Color online) Comparison of suppression between D, 
B mesons and non-prompt Jjijj in 2.76 TeV Pb-Pb collisions 


of the suppression of D meson and non-prompt J/ip si¬ 
multaneously. 


IV. EVOLUTION OF HEAVY MESONS IN A 
HADRON GAS 

As has been discussed in the previous section, at the 
critical temperature Tc, both the QGP fireball and heavy 
quarks hadronize into color neutral bound states. We can 
obtain soft hadrons from the bulk matter via the Cooper- 
Frye formalism and obtain heavy hadrons through our 
hybrid hadronization model. Subsequently, the produced 
hadrons from each event are subject to hadronic rescat- 




FIG. 9: (Color online) Effects of hadronic interactions on the 
D meson (a) Raa and (b) V2 in 2.76 TeV Pb-Pb collisions. 


tering which is modeled through UrQMD [5l|, . 


Unlike the Langevin equation that only requires a sin¬ 
gle transport coefficient, the UrQMD model requires the 
microscopic cross sections of hadronic scatterings as cru¬ 
cial inputs. To simulate the rescatterings of D mesons 
inside a hadron gas, we introduce into the UrQMD frame¬ 
work the scattering cross sections for charm mesons with 
TT and p mesons as calculated in Refs. [toI - I^ which 
are based on a hadronic Lagrangian generated from local 
flavor SU(4) gauge symmetry. In this calculation, uncer¬ 
tainty remains in the choice of the cutoff parameter in the 
hadron form factors. We treat the variation in the cutoff 
as a systematic uncertainty in our following calculations 
of heavy meson observables. 


In Fig. 9(a)[ we investigate how D meson Raa is 
affected by the hadronic interactions. One observes 
that due to the additional energy loss experienced by 
D mesons inside the hadron gas, Raa for D mesons is 
further suppressed at large pt- Consequently, due to the 
conservation of the number of charmed hadrons, D me¬ 
son Raa is slightly enhanced at low px after the UrQMD 
evolution. As mentioned above, the error bands in our 
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FIG. 10: (Color online) Effects of hadronic interactions on the 
D meson (a) Raa and (b) V2 in 200 GeV Au-Au collisions. 


results characterize the uncertainties introduced by a fac¬ 
tor of 2 difference in the choice of the cutoff parameter 
in the hadron form factors when calculating the heavy 
meson scattering cross sections in Ref. 0. With our 
comprehensive framework that incorporates heavy flavor 
evolution in both QGP and hadronic phases, we provide 
a good description of the D meson suppression as ob¬ 
served in 2.76 TeV central Pb-Pb collisions. After the 
inclusion of the hadronic interactions, the spatial diffu¬ 
sion coefficient of heavy quarks extracted from high px 
Raa data is updated to 6/(27rr). 

The effect of the hadronic interactions on D meson V 2 
at the LHC energy is shown in Fig. 9(b) Due to addi¬ 
tional scatterings of D mesons in an anisotropic hadron 
gas, its V 2 is further enhanced by around 20%. In Fig. 
|9(b)[ we also present the difference between two hydrody¬ 
namic initial conditions. Since the KLN model provides 
a larger eccentricity of the initial entropy density profiles 
than the Glauber model, this may cause another 20% dif¬ 
ference in the collective flow of heavy mesons after their 
evolutions inside the QGP and the hadron gas. How¬ 
ever, after taking all effects into account, our calculation 
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FIG. 11: (Color online) The D meson suppression in different 
centralities at the RHIC experiment. 



still underestimates D meson V 2 compared to the ALICE 
data. Several studies has been carried out targeting this 
V 2 . puzzle. For instance, it has been suggested in Refs. 
[Ta [74| that by taking into account the temperature de¬ 
pendence of the transport coefficient (q/T^) and increas¬ 
ing the relative contribution of the medium modihcation 
to heavy flavor spectra around Tc, the anisotropy pa¬ 
rameter V 2 in the final state can be effectively enhanced. 
These effects will be investigated in detail in the future. 

In Fig. uni we provide our calculations of the D meson 
nuclear suppression factor and anisotropic flow param¬ 
eter for Au-Au collisions at the RHIC energy. Similar 
to the LHC scenario, the hadronic interactions simulated 
with the UrQMD model slightly suppress D meson Raa 
at large px and enhance the anisotropy parameter V 2 - 
Our numerical results are consistent with the experimen¬ 
tal data from the STAR Collaboration. 

In Fig. El we present D meson Raa for different cen¬ 
trality classes as measured in the RHIC experiments. In 
Fig. m we show the integrated Raa of D mesons over 
given Px regions as a function of centrality characterized 
by the participant numbers. One can see that as moving 
from more central to more peripheral collisions, D me¬ 
son Raa increases due to a smaller geometric size and 
a shorter lifetime of the hot and dense nuclear matter 
created in more peripheral collisions. Our calculations 
are consistent with all the available data from the RHIC 
and a prediction for the participant number dependence 
of the D meson Raa is also provided for a smaller px re¬ 
gion. We note that the difference of the integrated Raa 
between the 0 < px < 3 GeV regime and 3 < px < 8 GeV 
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FIG. 12: (Color online) The participant number dependence 
of the D meson Raa at the RHIC. 


is expected to depend on the centrality of collisions due 
to a combination of heavy flavor energy loss and the co¬ 
alescence mechanism in heavy meson production. 


V. SUMMARY AND OUTLOOK 

In this work, we have established a comprehensive 
framework to describe the full evolution history of heavy 
quarks together with the evolution of the fireball in rel¬ 
ativistic heavy-ion collisions, including their initial pro¬ 
duction, energy loss in a QGP medium, hadronization 
an the subsequent interactions of heavy mesons with the 
hadron gas. At the beginning, the entropy density of 
the bulk matter produced in nuclear collisions is initial¬ 
ized with either the MC-Glauber or the MG-KLN model; 
heavy quarks are initialized with the MG-Glauber model 
for their position space distribution and a pQGD cal¬ 
culation for their momentum space distribution. Dur¬ 
ing the QGP stage, the bulk matter evolves according 
to a (2-|-l)-dimensional viscous hydrodynamic model, 
while the heavy quark transport inside this medium is 
described by our modified Langevin equation incorpo¬ 
rating both quasi-elastic scattering and medium-induced 
gluon radiation processes. At the critical temperature Tc, 
the QGP matter is converted into soft hadrons accord¬ 
ing to the Cooper-Frye formalism, and heavy quarks on 
the other hand hadronize based on the hybrid model of 
fragmentation and coalescence model we develop. In the 
last stage, both soft and heavy hadrons are fed into the 
UrQMD model for the simulation of hadronic scatterings 
until all interactions cease. Our numerical framework is 
designed in such a way that each evolution stage can be 
easily replaced by another model, e.g., a different hy¬ 
drodynamic background, a different heavy quark trans¬ 
port model or a different hadronization process, therefore 
a systematic comparison between different theoretic for¬ 
malisms can be conveniently implemented in the future. 

With our current approach, we have shown that while 


the collisional energy loss dominates the low px region of 
heavy quark transport inside the QGP, the contribution 
from the medium-induced gluon radiation is significant 
at high pt- During the hadronization process, the frag¬ 
mentation mechanism dominates the high px regime, but 
the introduction of the heavy-light quark coalescence sig¬ 
nificantly enhances the production of heavy meson at in¬ 
termediate px in nucleus-nucleus collisions. In addition, 
the hadronic interactions after the QGP decays further 
suppresses the heavy meson Raa and enhances its ellip¬ 
tic flow V 2 ■ In this work, the mass dependence of heavy 
flavor evolution is also investigated. It has been found 
that due to the larger mass of bottom quarks compared 
to charm quarks, the former lose less energy. The effect 
of such a mass hierarchy on the final heavy meson spec¬ 
tra fades away around px = 20 GeV for the collisional 
energy loss, but still remains up to 40 GeV for the radia¬ 
tive energy loss. Also due to the larger masses of bottom 
quarks, the coalescence dominates over a wider px re¬ 
gion for the hadronization as compared to charm quarks. 
Within our framework, we have provided numerical re¬ 
sults of the heavy meson suppression and anisotropic flow 
coefficients that are consistent with most data from both 
the LHG and the RHIG experiments. The spatial dif¬ 
fusion coefficient D of heavy quark extracted from our 
model is between 5/(27rT) and 6/{2ttT), depending on 
whether the hadronic interaction is included or not in 
the calculation. These numbers may be translated into 
qa/T^ around 9.4 ^ 11.3 for a gluon jet, or qp/T^ around 
4.2 ~ 5.0 for a quark jet, which are consistent with the 
value obtained by the JET Gollaboration by comparing 
various light flavor jet energy loss formalisms to the ex¬ 
perimental data [^ . 

Our study constitutes an important contribution to¬ 
wards a more quantitative and accurate understanding 
of the full evolution of heavy flavors produced in rela¬ 
tivistic nuclear collisions. Several further improvements 
await our future effort. For instance, the nature of heavy 
quark dynamics in the pre-equilibrium stage of heavy- 
ion collisions [zllzi is still not clear at this moment, 
which may affect the final state hadron spectra. It has 
been suggested that heavy quarks produced by the gluon 
splitting process may experience different medium mod¬ 
ification pattern compared to those directly produced 
through the hard scatterings As discussed earlier, 

by increasing the relative contribution of energy loss near 
Tc, the heavy quark V 2 may be increased a lot without 
affecting its overall suppression [tI, this might be 
helpful to explain the large V 2 puzzle. These aspects will 
be explored in our future work. 
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